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Glossary of Symbols 


c 

D 

e 

E 

F 

h 

K 

m 

M 

n 

r 

R 

t 

T 

V 

X 

V 

P 

o 


specific heat 

2 -1 

diffusion coefficient, cm sec 
vapor pressure, dynes cm ^ 
equivalent grey body emissivity 
correction factor (Fuchs, 1959) 
enthalpy of solution per unit mass 
thermal conductivity of air 
nid.s s / qth 

molecular weight, gm per mole 

number density of particles or droplets, cm 

radius of droplet, cm 

specific gas constant, ergs gm ^ 

time , sec 

temperature, ®K 

ventilation coefficient (Squires, 1952) 

mixing ratio, gm per gm of dry air 

molality of HCl^q 

molecular viscosity of air 
-3 

density, gm cm 

surface energy, ergs cm ^ 

-2 -1 

Stef an-Boltzmann constant, ergs cm sec ®K 


Subscripts 


a 

ambient air 

s 

saturation 

i 

size category 

V 

vapor 

k 

heat 

w 

water 


r droplet or droplet surface 



I . Introduction 


In approaching the problem of the drift, dilution, 
removal and deposition of the SRM exhaust products, we face 
complexities that are in major part attributable to the 
very great range of the scales of atmospheric processes 
(Table 1.1) In modeling by numerical techniques this consi- 
deration is crucial, because a given numerical model can 
span only a limited domain to which the scale (and hence the 
associated "time constants” or "decay periods") must 
be restricted. 

Specifically, the exhaust plume problem starts with 
processes of molecular scale (10 ^cm) as particles form in 
the hot concentrated exhaust gases, and we are plunged 
immediately into problems of coagulation, diffusion and 
sedimentation of particles (10 ^cm to lO^^cm) and gases, 
including air and water vapor entrainment, in forming what 
we call the SGC. The meteorological problem, by mutual 
understanding, begins with the SGC, some 90 seconds after 
ignition: a ball of '\-250 m radius at about 500-700 m 

altitude, full of earth surface debris in addition to the 
HCl, AI 2 O 2 , and H 2 O specifically produced by the roclcet 
motor combustion. 

We are concerned about the disposal of this SGC, and 
we recognize at least 3 specific kinds of processes; (1) sed- 
imentation and dry deposition, i.e., the fallout of large 
-3 

(r > 10 cm) particles, plus the turbulent and diffusive 

-4 

impaction of smaller particles (r > 10 cm) on surface 

_3 

obstructions, and the direct sorption of gases (r ''' 10 cm) 


2 


by plants, etc; (2) scavenging by rainfall, which we shall 
expand upon later; and (3) diffusion, dilution and transport 
by atmospheric motions until background concentration levels 
have been reestablished. 

In particular, the present effort is especially 
focussed upon those removal processes that involve cloud 
formation and rainfall. We have launched two parallel 
efforts in attacking this problem, one at the small scale 
end of the spectrum, which is called "cloud microphysics", 
and the other at the mesoscale (Table 1.1), which in effect 
sets the stage within the atmospheric environment for the 
development of cloud and precipitation. 



Table 1.1 


Scales of Meteorological Systems 
and Processes. 


Synoptic Systems 
Mesoscale Systems 
Convective Showers 
Individual Clouds 
Air Turbulence 
Rain Drops 
Cloud Droplets 
Nuclei 

Mean Free Path 
H 2 O Molecule 


1 

'V'lO cm 

= 

lO^km 

7 

'\^10 cm 

= 

lO^km 

1 

'v-lO cm 

= 

10 )cm 

'\'10^cm 

= 

1 km 

H 

0 

- 

10°cm 

10 ^ 

- 

10 ^cm 

10“^ 

- 

10 ^cm 

10-6 

- 

10 ^cm 


'\i 

10 ^cm 

3 

X 

10 ^cm 
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II. Cloud Microphysics 

Cloud microphysics encompasses the processes of cloud 
droplet nucleation, growth and coalescence, which; serve 
to incorporate the constituents of the local atmosphere into 
the cloud elements, and to convert cloud to precipitation. 


A. The Implicit Model 


The diffusive growth of a droplet of pure substance 
from its vapor phase may be expressed by the modified 
Maxwell diffusion equation: 


dm _ 4rrrDF^V 
dt R 




( 1 ) 


In this expression 

is the correction factor due to Fuchs (1959) 
which resolves the continuum vapor flux with 
that given by kinetic theory at and near the 
surface of small droplets; and 
V is the ventilation coefficient of Squires (1952) 
adapted from Frossling (1938). 

All other symbols follow their customary meteoro- 
logical usage, and are defined in the glossary, p. 4. 
To complete the basic model of droplet growth, the 
heat transfer from the growing droplet must be considered and 
balanced against the latent energy released; 

L ^ = 4irrKFj^V(T^-T^) (2) 
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Mason (1957) devised a procedure involving a series of 
approximations aiid assumptions, whereby (1) and (2) are 
combined into a single growth equation that depends only 
upon the saturation ratio^ S = e/e_ (T ) , and the air 
temperature , T : 

cl 

^ _ s-1 
dt a+b 

where a = L^p /(KRT^); b = p RT /e (T )D. 

The growth equation, (3) has been widely used for 
the computation of cloud droplet growth upon particles 
composed wholly or partially of soluble inorganic salts. 
Because the droplet temperature, T^, is not explicitly 
required to solve (3) , this system is referred to as the 
implicit model. 


(3) 


B. The Explicit Model 

In its simplest conception, the HCI-H 2 O system with 
which we must deal in the case of the solid rocket motor 
(SRM) exhaust products, involves two volatile components 
and their interactions. Clearly, therefore, the microphysical 
problem must be reexamined. 

1. The Diffusion Equations 

When HCl and H 2 O vapors are present together, they 
jointly determine the equilibrium vapor pressure of each 
component over an HCl^q droplet, and Mason's principal 
simplifying steps are inapplicable. It is assumed that the 
two vapors may be treated as if they diffuse independently, 
so that the change of mass of a solution droplet may be 
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( 4 ) 


expressed as 

diti^ = dnij^ + dm2 

where the subscripts indicate the respective components, 
i.e., H 2 O and HCl. Then, analogously with equation (1) 


and 


dmj^ 

dt 


4-rrrD,F V 


R 


1 




dm 2 

dt“ 


4nrD_F V 
2 V 2 




(5) 


(6) 


2- The Energy Balance 

The energy balance of the droplet may be expressed by 




where the Q's are defined as follows: 

3T 


= m c 3^7— 
T r r dt 


internal energy increase 


(7) 

(7a) 


dm^^/dt; 


latent energy release (7b) 




= -4irrKFj(,V(T^-T^) 


conductive heat transfer (7c) 
from drop 


Q 


R 


-16Trr^ (Tj.~T^) 


radiative transfer from 
drop 


(7d) 


^ . 2,20 , 1 ^"‘r 

Qr = t-r at- 


3p 'mj^ dt 


W 


1 '^"'2 

- -= —) > 

m 2 dt 


-y 


ir,JL ^ 

3y dt 


^ dm2^ ‘^'*’r,3o 2o 1 

” ' dt ^3T “ 3 


aPr 


m 2 dt 


)] 


surface increase energy (7e) 
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mixing 


(7f) 




= 


dm 

- c (T -T 

r r a dt 

- 2 2 
2r m^g 

9iTZ 


frictional conversion (7g) 


, ,, dm 

y dh r 

m dy dt 
r 


heat of dilution (7h) 


The latent energy term (7b) accounts for phase changes of 

both H 2 O and HCl in the present context. Although L 2 is 

about 17 kcal mol ^ as against an value of about 10.7 

kcal mol at most stages of growth *^^1 ^ *^^2 and the 

dt dt ' 

heat requirement/contribution of HCl phase changes may be 

comparable to that for H 2 O. The radiative transfer term 

(7d) depends upon the value of the equivalent grey body 

emissivity, E, which is not precisely known. We estimate 

that it lies within the range of 0.7 to 0.9 or so for the 

infrared region of concern. The surface energy term (7e) , 

because of dependency of both solution density and surface 

tension upon molality and temperature, is quite complicated. 

The mixing term (7f) accounts for condensation of vapors 

at air temperature, T , and subsequent adjustment of the 

liquid increment to the droplet temperature, T^. For an 

evaporating droplet, is assumed to be zero. The friction 

term (7g) accounts for the conversion of potential energy 

to heat by dissipation at terminal fall speed. It is 

obviously very small (see below) for cloud droplets. The 

dilution term (7h) accounts for the requirement/contribution 

of heat as the HCl _ solution changes in concentration. 

aq 
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Enthalpy of dilution is tabulated in standard chemical 
references (e.g. , • Weast, 1971). 

3. Relative Values. of the Energy Terms 

The system (4), (5), (6), (7) may be regarded as 
being applicable to any number of volatile components, as 
long as one allows that additional energy terms for interactive 
heat sources and sinks may be added to (7) as written above. 

In the present work we are directly concerned with the 
H 2 O-HCI system involving also the solid AI 2 O 2 exhaust 
products of the SRM. We shall therefore de-emphasize the 
broader generalization of the present scheme and focus 
upon its application to the SRM exhaust problem. 

To simplify (7) for the immediate purpose, it is useful 
to make an order of magnitude comparison among the respective 
energy balance terms. Obviously Q and Q are the largest 
terms. They have opposing effects upon the droplet temper- 
ature, T^, and are effectively differenced in the process 
of computing T^. In fact, the fundamental assumption that 
permits elimination of T^ from the implicit model growth 
equation is that Q, and Q are identically equal. We note, 
however, that this is only approximately true, and that, 
indeed, the actual growth rate is necessarily dependent 
upon the residual between Q,. and Q . Terms in the energy 

Ij 

balance equation that are of the same magnitude as this 
residual must therefore be considered non-negligible for 
the explicit model. 

Utilizing output from implicit model calculations of 
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the growth of cloud droplets on a population of (NH^) 2 S 0 ^ 
nuclei in a 1 m sec ^ updraft 20 m above cloud base, 
estimates of the respective terms of (7) have been computed 
and are presented in Table 2 as multiples of Q^, the internal 
energy (heat storage) term. It is clearly shown for this 
case that and are 3 to 5 orders of magnitude larger 
than for droplets up to 10 ym radius. This ratio 
diminishes as drop size increases, and will be much smaller 
for 100 to 1000 ym-radius raindrops. Terms within one order 
of magnitude of are Qj^ and Q^, each of which diminishes 

with increasing droplet size. Inasmuch as is an increasing 
function of droplet size, it contributes to the diminution 
of the above ratios as droplets become larger. The absolute 
values of are therefore given at the bottom of the table 
in units of erg sec The frictional term increases with 

droplet size inasmuch as it depends mainly upon droplet 
terminal fall speed. Obviously the values for cloud drop- 
lets are so small (11 to 15 orders of magnitude less than 
as to be totally negligible for modelling purposes. 

We note however that fall speeds for these droplets are of 
order 10^ less than those of raindrops, for which therefore 
the frictional term should be about 10^ larger. The dilution 
term, Qj^, is also quite negligible for (NH^) 2 SO^ in cloud 
droplets, but may become marginally significant for raindrops. 

The dilution heat for HCl is much larger. For 4.0 molal 
aquaeous HCl droplets, Q^/Q^ takes values of 370, 29 and 11 
respectively for droplets of radius 1, 5 and 10 ym growing 
at the rates derived from the above implicit model calculations. 

Thus, for the HCL/H 2 O system, the heat of dilution is substantial. 
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Table 2'.'2 


: Relative Values of the Energv 

Terms for { (NH^) 2S0^} . * 



r * 1 ym 

r = 5 ym 

r = 10 ym 

Ql/Qt 

5.913E4 

4.461E3 

1.719E3 

Qr/Qt 

5.911E4 

4.456E3 

1.716E3 

Qr/Qt 

16.0 

4.92 

3.55 

Qc/Qt 

3.34 

5.16E-2 

1.02E-2 

Qm/Qt 

2.62 

4.40E-1 

2.53E-1 

Qp/Q^ 

1.59E-15 

2.58E-14 

9.13E-12 

Qq/Q^ 

1.07E-6 

3.82E-4 

5.9E-3 

For 4 . 0 

molal HCl ; 

aq 



Qjj/Qt 

370 

29 

11 


♦Nuclei vertical rise rate of 1 m sec 



20 m above cloud base. 
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The contribution of Q- deserves specific mention. 

2 

The molar latent heat, L 2 , released as HCl makes the vapor- 
liquid transition of the order of 1.7 x Lj^. Because this 
heat warms the drop, it must operate to slow droplet 
growth of the relatively pure H 2 O droplets in a mixed 
population having a uniform environment, as they acquire 
HCl from the vapor more rapidly than their neighbors. 
However, as the SGC progresses in time, the HCl partial 
pressure is expected to decrease rapidly in response not 
only to the solvation in droplets, but also because of 
turbulent dilution of the cloud by environmental air. It is 
expected that, at some point, the dilution will reverse the 
HCl vapor pressure gradient, whereupon HCl will begin to 
leave the liquid phase, and L 2 will be ‘extracted from, 
rather than released to, the aquaeous HCl droplets. The 
cooling effect of this process must then contribute to an 
increased condensation of H 2 O from the air, and may be 
an important mechanism for accelerating the onset of rain. 
That this capability is particularly a property of the most 
strongly acidic droplets is also of particular significance 
in determining the deposition pattern of aquaeous HCl. 

An additional attempt to adapt the implicit model 
output to the study of the HCI/H 2 O system was made. In 
this case the values after 10 sec of upward displacement 
from cloud base at 1 m sec ^ speed were adapted from the 
(NH^)2S0^ calculations. The model values for dm/dt, m^, 
and dT^/dt for (NH^)2S0^ nuclei, (T^-T^) , and y are 
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r = 0.1 ym 

1.0 ym 

10 ym 

dm/ dt 

9. 4E-16 

1.66E-12 

9.8E-11 

m 

4.42E-15 

3. 83E-12 

4. 96E-9 

w 

dT^/dt 

8.6E-3 

6. OE-3 

6. OE-3 

T -T 

2. llE-4 

1.46E-2 

7. llE-2 

r a 

V 

1.55E-1 

4.4E-3 

9.8E-2 

at 10 m above 

cloud base. 

Using these 

values as if 


were directly applicable to condensation, the 

figures in Table 3 are computed. The droplet of r=l ym 

grows more rapidly relative to its mass than the 0.1 ym and 

10 ym droplet. It is also the most dilute of the droplets 

presented (see above) . The dilution term is similar in 

relative magnitude to its value at 20 -m above cloud base 

(above). For other HCl concentrations, multiply 
2 

Qd/Qt by (y/4) . clearly the relative magnitudes of the 
respective energy terms depend upon droplet size, rate of 
growth and chemical species; and different stages in the 
processes of condensation, evaporation, coalescence and 
precipitation may permit neglect of or may require retention 
of different terms in the energy budget. Only Q_ appears 
to be totally negligible, and this allows simplification of 
(7) to 

Qt = °L Qr - % Qm ^ Qd 

By means of the definitive expressions (7a-h) applied to (8) , 
the droplet temperature elevation may be written 
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Table 2v3: Estimates of the Energy Terms in 

Relation to Q_ for the HCl/H^O System 
(see text) af^y = 4. ^ 



r - 0.1 ym 

r = 1.0 ym 

r = 10 ym 

Ql/Qt 

1.46E4 

4.28E4 

1.95E3 

Qr/Qt(e=1) 

0.92 

9.7 

3.7 

Q^/Qt 

8.9 

2.6 

0.012 

Qjj/Qt 

0.006 

1.05 

0.24 

Qp/Q^ 

8E-15 

1.2E-12 

1.2E-10 

Qp/Q^ 

94 

276 

13 
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(9) 


(T^-T^) = [- 




[4Trr(KFj^V+4raj^ETp + C^. dmj./dtl 


The diffusion equations (5) , (6) , and the droplet 

temperature elevation equation (9) represent the micro- 
physical processes that require evaluation within each time 
step for each size category of the cloud droplet spectrum. 
The size spectrum is usually well-enough expressed by 25 
to 50 size categories. 

In addition, for each time step, the mass conservation 
accounts are kept. Assuming no dilution of the air parcel 
by turbulent mixing with environmental air 



i i^dt V 


( 10 ) 


_ E .‘^"' 2 . 
dt i "^i ^dt ^ i 


( 11 ) 


The enivornmental temperature and pressure changes from one 
time step to another are given by 


dt 


/ . dw. 

'"Pa^^ dt> 


(12) 


and 


dT 

C 

dF 


dx, dx„ , , 

^^1 dt~ ^2 dt~ ^^'*'dt^ dl'^^'^^l'*'^2 ^ ^ 

'^p ^1^P3^ ^2% i^^^r^r^i^i 


(13) 


where the summation in the denominator represents the heat 
capacity of the liquid carried with the parcel. Equations 
(5), (6) and (9) for each droplet size category, plus 
equations (10) , (11) , (12) and (13) constitute the explicit 
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model for two volatile components. 


4. Solution of the Explicit Model for Two Volatile Components 
a. 


Initialization 
,( 0 ) 


The values of T^'" , , , n (r^) , and w^®^ 


are specified for time t^^^ = 0 by the statement of the 

problem. The values of T^®^ may be determined for each 

droplet size category, if m. and mf®^ are known, as 

^1 ^2 

follows : 

1. Assume an approximate value of T 


( 0 ) 


(0)^ _ ^(0) 
i 




where { } designates an estimated value. 
2. Calculate the corresponding values of 


dm(0> dm^O) 

] and [ 


] for each size category 


'dt " "dt 

using (5) and (6) . Here the square brackets [ ] 
designate a calculated value. 

Calculate [T^®^] by means of (9) using the results 
of step 2. 


Compare the estimated against the computed 

[T^^^], relax the difference and recompute until 
^ i 


,( 0 ) 


{"^r. } - [T^°h < e, 0(10 


-10 


The values of mf®^ and m.^®^ are derived from measure- 
il X2 

ments of the SGC. For example, one may assume following 
Squires (1976) that each AI 2 O 2 particle adsorbs a monomolecular 
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layer of HCl, leaving the HCl residue as vapor to be equili- 
brated with the condensing H 2 O coating on each particle. 

The role of the adsorbed HCl in determining molality 
of the aquaeous phase may be modeled in various ways 
(i.e., assuming some proportion of solution and some 
residual adsorption bonding) in the absence of hard data. 

We depend upon our chemical advisers for aid in designing 
this part of the model. 


b. Subsequent Time Steps 

For subsequent time steps, = t^^^ + At^^^ + 

At^^^ + ... + At^^ we calculate dT^^^dt using equation (13). 

For values of dxj^/dt and dx 2 /dt we use those determined 
from the prior time step. To compute dT^^V^t we set 


dx<-l> dx<-") 


= 0 


dt dt 

For the droplet growth computations we define the 


parameter 


k-2) 


(T -T ) (T -T ) 

r a r a 


At 


and 


(k-2) 
(k) 


use it to predict { (T -T ) ' ' } by means of 

^ cl 


(14) 


{ (T -T ) ^^^ } = (T -T ) At^^“^^ 

r a r a 

For k = 1, we define = 0. This is equivalent to 

setting 

{ (T -T^) } = (T^-T ) 

which is a reasonable first estimate. The steps 2, 3, and 

4 (see above: a. Initialization) are repeated to obtain the 

final values of T^^^ and dm, o^^^/dt. 

r 1, 2 ' 
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5. Numerical Procedure 

The predictor-corrector method of Hamming (Ralston 
and Wilf, 1960) is used for general integration of the 
explicit model at present. This technique does not "self- 
start" ; so a fourth order Kutta-Simpson procedure is used 
through the first two time-steps. Although it is quite 
likely that more efficient procedures may be applied, our 
present orientation is toward generating a working model. 

We expect to take up the computational efficiency question 
separately at a later date. 

6. Verification 

Examination of the model derivations indicates that (a) 
the assumptions and approximations used to obtain the implicit 
droplet growth equation should lead to slightly lower than 
true growth rates; (b) the explicit equations should produce 
faster growth rates than the implicit equation gives, but 
their relation to the "true" growth rates is not easily 
judged. 

Inasmuch as both systems can be applied to water 
condensation upon nuclei composed of inorganic salt, we 
have modeled the early stages of cloud droplet growth 
by each model using identical nuclei spectra and external 
conditions. The results for ammonium sulfate nuclei are 
shown in Fig. 2.1 for a steady updraft of 1 m sec ^ to a 
height of 100 m above cloud base. Here the explicit model 
gives slightly larger droplets than the implicit model. 


18 



DROPLETS PER CM^ PER MICRON RADIUS INTERVAL 



RADIUS, /xm 


Fig. 2.1: Droplets per cm per Micron Radius Interval- 

100 meters above cloud base, steady 
updraft. 
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Figure '2.2 shows results from both models for a simulated 
parcel trajectory that rose 20 m above cloud base at 
1 m sec then returned to cloud base at 1 m sec ^ and rose 
again at this same rate to 100 m above cloud base. The 
nuclei again are ammonium sulfate, and both models give a 
bimodal distribution after 140 sec. The greater differences 
of size between the explicit and implicit model results is 
attributed partly to the longer time interval and partly to 
the failure of droplets smaller than the fourth size category 
to reactivate after the 20m-amplitude vertical displacement 
cycle in the explicit model calculations. During the second 
upward displacement, the saturation ratio does not achieve 
as large a value as during the first because the presence 
of the larger sized droplets provides' an increased vapor 
sink as the second uplift stage begins. 

We conclude that the explicit computations are adequately 
validated by the above for (NH^)^SO^ nuclei. Additional 
tests are planned for salts having larger values of the 
heat of solution. NaOH, KOH and KCIO^ appear to be likely 
candidates for these computations. 
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III. Meso-Scale Model 


A. Introduction 

A meso-scale meteorological model is being developed 
to study the diurnal and seasonal variations of the meso- 
scale weather patterns in the lower atmosphere over the 
Florida Peninsula. The basic structure of the model follows 
the work of Pielke (1974) . We have introduced modifications 
for the purpose of removing some of the approximations 
of the original model. In particular parameterized effects 
of cloud and precipitation processes are introduced, and the 
boundary-layer processes are treated in more detail than 
in Pielke 's three-dimensional model. Nonetheless, these 
processes are still highly parameterized, because they cannot 
be resolved explicitly in the meso-scale. The detailed 
formulations of this model are presented in the following 
sections. 

The results from this model will be used as continuous 
inputs to a sub-meso-scale model for the area of Cape 
Canaveral to predict the processes and patterns of SRM 
exhaust product removal after Space Shuttle launches. 

B. Basic Equations 

Several basic assumptions are made in deriving the 
system of equations that express the atmospheric processes. 
They are 

(1) that the ideal gas law is valid for moist air, i.e. 
p = pRT^ 
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( 2 ) 


and 


that the shallow atmosphere has constaht density, 

i.e. 


dp ^ 
dt 


C3) that a hydrostatic balance is held for the meso- 
scale motions, i.e. 


= - 
3z 


pg» 


where p is the pressure, p the density, the virtual 
temperature, R the specific gas constant for the air and 
g the gravitational acceleration. 

After performing the scale analysis, the basic 
equations in Cartesian coordinates for the meso-scale 
motion are. 


du 

dt 


- If + 


+ F 


u 


( 2 . 1 ) 


dv 

dt 


- i - fu 
P ciy 


+ F 


( 2 . 2 ) 


1 ^ 
p 3z 


pg 


(2.3) 


9u . 3v . 3w 


0 


(2.4) 


p = pRT^ 


(2.5) 


^ - -1-fS + P 1 

dt CpT'^0 ^ ^6' 


^ = S + F 

dt Sq + Fq 


( 2 . 6 ) 

(2.7) 
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where 





3 


+ V 


a 

ay 


+ w 



(i) 


F 


Ci) 


^ ax 


a(ii 

ax 


) + 


a 


IK 


H 


3(i) x 

ay ‘ 


+ -2-(k 

3 Z ' z 


(i) 


a(i) 

iz 


)] 


f = 2il sin 
f* = 2J2 cos <j>. 


u, V, and w are the x-, y-, and z- components of velocity 
respectively, 0 the potential temperature, q the specific 
humidity, n the earth rotation rate, and 4> the latitude. 

K„ and are the horizontal and vertical eddy exchange 

n Z 

coefficients. K„ has the following form (Leith, 1969) ; 
n 

Kr = C^(AX) (Ay){ (|^ + 1^)^ + l/2t(|J)^ + (|~)^]}^/^ 

where is a constant. Ax and Ay are the x- and y- direction 
grid spacing. is defined below (section 6) in the 

discussion of the boundary- layer parameterization. S„ 

W 

and Sq are the source functions respectively for the potential 
temperature and the specific humidity. The details can be 
found below in section 5 on cloud and precipitation pareuneter- 
ization. 

Now we introduce two dependent variables. First, the 
non-dimenaional pressure (tt) has the following definition, 

R/c 

ir = Cp(p/p*) P (2.8) 


where Cp is the specific heat capacity at constant pressure, 
and p* (=1000 mb) the reference pressure. Secondly the 
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defined 



aasic equations may be 


( 2 . 10 ) 

( 2 . 11 ) 

( 2 . 12 ) 

( 2 . 13 ) 

( 2 . 14 ) 

( 2 . 15 ) 


( 2 . 16 ) 



C. Boundary Conditions; 


To eliminate reflection from the boundaries, and to 
have smooth continuation of the flow through the boundaries, 
the following boundary conditions are assumed: 


Cl) Lateral boundaries: 

= ^ = M - n 

cix 3x '3x 3x ~ ^ 

3u _ 3v ^ ^ 

~ 3y 5y ° 


3x 


3y 


at X s= 0-, and L^; 
at y = 0 f and L 2 , 


where and L 2 are the lengths of the horizontal model 
region in x- and y- direction respectively. 

( 2) Bottom boundary : 

The bottom surface is flat, and u = v = w =0 . 

The virtual potential temperature (8^) and the specific 
humidity (q) are determined by time-dependent functions which 
simulate the diurnal variation and also imply solar and 
nocturnal radiations. These functions are expressed by 
truncated Fourier series. They operate to perturb the 
initial balanced state and indirectly to generate the meso- 
scale motion. 


D. Initial Conditions 

The following data are needed to initiate the model : 

(1) a three-dimensional moisture field, 

(2) a three-dimensional temperature field, and 

(3) a two-dimensional surface pressure field. 

The specific humidity, q, is derived from dew-point, temperature 
and pressure data. The virtual potential temperature and 
the non-dimensional pressure are calculated by (2.8), (2.9), 
and (2.12). 
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The initial horizontal velocity field is determined on 
the basis of two assumptions; first, that the lowest 5 
levels in the model are assumed to have the Ekman solutions 
(Holton, 1972) , and, second, that geostrophic balance pre- 
vails at the higher levels. The initial vertical velocity 
is computed by (2.13). 

The initial height of the top free surface is given by 
^ = ^Psfc ■ Ph^/P9. 

where is the surface pressure, and p^^ = 500 mb. 


B. Cloud and Precipitation Parameterization. 

The effect of latent heat release from clouds and 
precipitation upon the meso-scale motion is parameterized 
into two major parts. This allows the source terms in (2.14) 
and (2.15) to be expressed by 


S 

S 

q 



+ S, 


+ S 


and 


(5.1) 

(5.2) 


where subscripts 1 and 2 denote stable and convective 

precipitation respectively. 

(i) Stable precipitation 

By assuming that the stable precipitation 

process is pseudo-adiabatic, S„ and S may be written in 

»i q^ 

the form; 
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} (5.3) 


^ ^-L dq^/dt, when q>qg 
01 0 , when w£.0 


and w>0. 


S =r dq /dt, when q>q and w>0,-, 
q^ I s s t 

0 , when qfcjg or w<^0 


(5.4) 


where q^ is the saturation specific humidity, and L the 
latent heat of condensation. 

The rate of change of q^ is given by 




c R T 
P V 


T-eL 
^ + q L 


q w 


(5.5) 


where e = 0.622, and is the gas constant for water vapor. 

The detailed derivation can be found in Haltiner (1970) 
with minor transformation to adapt to the present coordinate 
system. 

(2) Convective precipitation. 

The parameterization of convective precipitation for 
larger scales has been organized by Fraedrich (1973) using 
the concepts of compensating subsidence (Arakawa, 1971) and 
of lateral mixing (Kuo, 1965) . The source terms for 6 and q are 


S 


0 


2 


mCp(r^-r) + (l-a)pCpa(T^-'T)/At, 


(5.6) 
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and 



where m(z) is the vertical mass flux, the dry adiabatic 
process rate, F the environmental lapse rate, a the fraction 
of a horizontal grid area covered by convective cloud, and 
a the ratio of evaporation to condensation. The subscript c 
denotes saturated values on the moist adiabat through the 
cloud base and At is the time step of the numerical integra- 
tion scheme. 

The vertical mass flux m(z) can be determined by an 
entrainment function suggested by Fraedrich (1973) , 


1 3m(z) _ 0.4 o^ 

m(z) 3z ~ D 

c 

where is the cloud diameter, and it is linearly related 
to cloud depth (Zq) (Plank, 1969-; Belts, 1973), i.e. 


D 

c 


z 


0 




1 . 


(5.9) 


The fraction, a, of the grid area covered by convective 
cloud is the ratio. 


a 


C2+C3 


(5.10) 


Here is the total integrated horizontal moisture conver- 
gence in the layer between the cloud base (Zg) and the level 
(z^) where the moist adiabat through the cloud base intersects 
the environment sounding, plus the surface evaporation, that is 
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= AT { - j 2 ^P[| 3 f(qu) + |yCqv)]dz + } , (5.11) 


where 

^v'PVo'W 

Cjj is the drag coefficient, the wind speed at the level 
just above the surface, the specific humidity at the 

underlying surface, and -q^ the specific humidity at the 
level just above the surface. 

C 2 and are expressed by 


C = I 

2 I Zg p(q^-q)dz, and 

r . <= 

Jzg Pj^(T^-T)dz. 


(5.13) 


(5.14) 


The sum of C 2 and represents the total water vapor that 
can be condensed. 

Finally, the weighting factor a is the ratio of 
evaporation to condensation, and (1-a) is the ratio of rain- 
fall to condensation. 


F. Boundary-Layer Parameterization 

The atmosphere below 500 mb in this model is divided 
into three physical layers. They are (1) the surface 
boundary layer (SBL) , (2) the planetary boundary layer 
(PBL) and (3) the free atmosphere layer. The SBL is a 
constant flux layer attached directly to the surface. The 
PBL is immediately above the SBL, and the region above the 
PBL is the free atmosphere. The atmosphere above 500 mb is 
assumed to be stable. 
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which 


The vertical eddy exchange coefficient, 
appears in equations (2.1, 2.2, 2.6 and 2.7) may be assumed 
to be a function of z. Because of the different stability 
conditions between daytime and nighttime, it is difficult 
to formulate a universal function to represent a full diurnal 
cycle. For daytime, the formulae in O'Brien (1970) are 
adopted; for nighttime, a functional expression is being 
developed. 

A summary of the vertical profile of in the 

daytime is given as follows. 



+ [ (H*-z)^/(H*-h*)^H ) + 

(z-h*) 2 (K^^^-K|^^)/(H*-h*)]}, forh*<z<H* 

{ (i) 

K£ ' , for H*<z 
(z/h*;K 2 ^^, ror z<h* 


where K^^^and are the values of at the 

altitude at the top of the PBL, and at h* , the altitude 

at the top of the SBL, respectively, 

(i) 


K 


(i). 


3k: 


(- 


9z 




( 6 . 2 ) 


and superscript i can be momentum, potential temperature, 

or specific humidity. is a constant in the free atmosphere. 

A Bermite interpolating polynomial suggested by O'Brien (1970) 

determines the distribution of K^^^in the PBL. In the SBL, 

xi^^is a linear function of z. 
z 


32 



H* is given by a time-dependent, empirical equation 


introduced by Deardorff C1974) , 

3 H* 3 H* 3 H, 

= - tu 


3t 


-^ + V -y^) + Wt + 
l.BCw* + 1. lU* - 3.3U*fH*) 


(6.3) 


where 


h; 
g — 
0 

V 


+ gw; + 7.2U^ 


H. 


w? = aw* + (l-a)w, 

X * 


(6.4) 


(- u*e*H*) if e*<o; 

W* ={ (6.5) 

0 , if 6 *> 0 , 

30V 

and (— ^)h is the virtual potential temperature gradient 
at H*, and U* and 6* are defined below. 

A relation between H* and h* is given by Pielke (1974) , 
and Pielke and Mahrer (1975) based on the work of Blackadar 
and Tennekes (1968) . That is 

h* = 0.04 H* (6.6) 

is determined by the micro-meteorological processes 
in the SBL. It is evaluated by means of the empirical 
expressions derived by Businger (1973) , on the basis of 
relatively recent field data (Businger, et.al. 1971). 
has the following forms. 

= kU*h*/(|.j^ (6.7) 
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Kj®’ - - kU.h,/*g 


( 6 . 8 ) 


where k is the von Karmcin constant. The dimensionless wind 


shear 


_ kz 3U 
'f’ M U* 32 ' 


(6.9) 


where 


n ^ 2 ^ 2,1/2 

U = (u + V ) ' 


and where the dimensionless virtual potential temperature 


gradient 


* - 2 
‘’’h 0* 3z 


(6.10) 


have been fitted to real observational data by Businger 
et.al. (1971). The empirical functions for them are. 


4>m { 


(1-155)“^/'^, C<0, 
1+4. 75, 5>70, 


( 6 . 11 ) 


0.74 (1-95) , 5<0, 

0.74 + 4 . 75 , 5>0, 


( 6 . 12 ) 


where 


5 = z/L, 

L = ®yg(uj/kz0*) , 


(6.13) 

(6.14) 


is the friction velocity, and 0* the scaling potential 
temperature. After integrating (6.9) and (6.10) with 
respect to z, we have 


U* = kU\ ^*/ [in {z/z^) -ip 
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and 


(6.16) 


0* = [0^ Ih* ~ ®v )ln(z/z^)-0.74ij;2^ ' 

2 £n[{l-*'(j)~^)/2] + Znia+^~^)/2]~2 tan"^((j)^^) + 
= { ir/2, 5<^0, 

- 4 . 75 , 5 > 0 , 


£n [ d+0.74(j)„^)/2] , 5<0 

^2 = ^ 

0.74 + 4.75 5>0 


(6.18) 


u|j^* is the wind velocity at h* , 9^|j^ and 9^i gf^, are 6^ 
at h* and at surface respectively, z^ is the roughness length 
which is different between water and land surfaces. Over 
water , 

a^ = 0.032 U*/g (6.19) 


suggested by Clarke (1970) is used, whereas over land it is 
assumed constant (4 cm in this model) . 

In order to close the system described above, one must 
know 5> The formulation of 5 involves the Richardson number 
Ri, which is defined by. 
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Ri = i 


3z 


6 (^)2 

V z 


(6.20) 


in the SBL, where 0 is the average virtual potential tempera- 

90V "" 

ture, - the virtual potential temperature gradient, and 
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3U 


-r— the wind shear. With the aid of (6.9), (6.10), (6.13), 

O Z 

and (6.14) we get an analytic relation. 
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Ri = 


(j)H 


U) 


M 


Substituting (6.11), and (6.12) 
algebraic equation for ^ , 


( 6 . 21 ) 

into (6.21) we have an 


where 


and 






+ 


d* = 


a* = -8.214 
b* = 0.5476 
c* = 9Ri^ 
d* = -Ri^ 


}for Ri<0 


a* = 0 

b* = 22.09 rI^4t7 
c* = 9.4Ri-0.74 
d* = Ri^ 


}for Ri>0 


( 6 . 22 ) 


The Richardson number in the present boundary-layer 
parameterization will be approximated by a quasi-local version 
of the Richardson number (rT) suggested by Portman (1974) , 

^^^"'l,2^®v2"®vl’ £nn*!/[e^ (U2~Uj^)2], (6.23) 

where 

" 1,2 - 

n* = h*/z^. 
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and subscripts 1, and 2 denote the values at h* and z^, 
respective ly . 

For the nighttime, the formulation for is being 

investigated, and will be presented later. 


S. Numerical Methods 

(1) Model structure 

A horizontal area (720 km x 480 km) is chosen to cover 
the Florida Peninsula in the center and the surrounding 
water (Figure 3^,i). The horizontal grid-space (Ax = Ay) is 
30 km. Study of grid spacing down to 15 km will be conducted 
to evaluate the effect of this parameter. The atmosphere 
below 500 mb will be divided vertically into ten levels 
with spacing (Table 3.1). The time increment (At) will depend 
on the computational criterion of the numerical scheme used. 

(2) Numerical schemes 

The differential equations shown in the previous 
sections will be approximated by the algebraic difference 
equations using finite-difference methods. The total 
derivative with respect to any dependent variable. A, 


dA _ 9A , 

dt at 




3A 

a 


(7.1) 
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720 km 
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Table 3 . 1 


Height of Fixed Vertical Level 


Level No. 

Height (km) 

11 

5.00 

10 

4.05 

9 

3.20 

8 

2,45 

7 

1.80 

6 

1.25 

5 

0,80 

4 

0.45 

3 

0.20 

2 

0.05 

1 

0.00 
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will be approximated by the forward-upstream scheme 
(Molenkamp, 1968) . It can be expressed as 


fdA. n 

'dt' . 

ipk 


1 ,,n+l 


' * ‘'’J jk + 0?Jk * ®?jk> ' <■'■2) 


where 




"?jk(^Jjk-A?-i,jk)/^^' “ijk> 0 


ijk 


, ,^-AV,^)/Ax, uV.vlO 


,n 


ijk' i+l,jk ijk 


n 


ijk- 




In (7.2) the subscripts i, j, and k represent the location 
of the grid points in the model. 


X = (i-l)Ax 
y = (j-l)Ay 
z - G(k) 

where G(k) is a function to space the uneven vertical levels. 
The superscript n shows the number of the time step. 


T = nAt. 

All other terms in the differential equations will be com- 
puted using the central-difference scheme. 
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The computational criterion for the forward-upstream 


scheme is 




<1 


where V* is the maximum speed of any component of velocity 
and A the corresponding grid spacing. 
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IV. Input Data on Equilibrium Vapor Pressures for HCl and 

H-0 over HCl 

Basic to the microphysical model calculations for the 
two volatile components, HCl and H 2 O, is the information 
on equilibrium concentrations of these vapors over HCl^^. 

For our purpose, it is important that the data reflecting 
these concentrations provide for the ranges of temperature 
and molality that are likely to be encountered in tropospheric 
cloud formation. It is a simple matter of circumstance 
that most chemical measurements available in the literature 
and handbooks have not been made with atmospheric applications 
in view. In this case, the best available data are those of 
Fritz and Fuget (1956) which give vapor pressures over bulk 
HCl^g for H 2 O vapor and HCl vapor equilibria in the range of 
temperature 273°K ^ T £ 323®K and for solution molality 
ranging from 0 to 15.88. 

Although initially, in dealing with the stabilized 
ground cloud (estimated altitude below 1 km) and the stabilized 
column cloud (altitude up to 3 km) over Florida, the temperatures 
to be encountered should not go much below 263 °K, we have the 
further consideration of ingestion into active convective 
cells in which much greater height and temperatures down to 
243°K or colder may be encountered. Thus it is necessary 
to extend the tabulated data of Fritz and Fuget (l.c.) 
beyond the limit of the measurements. 

An additional requirement upon these data for use in 
numerical modeling is that they form a smooth and continuous 
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functional surface throughout the domains of the data. In 
general, even with the best of experimental work, normal 
experimental error and accidental mistakes become incorporated 
into the results. It is a problem of statistical curve fitting 
and analysis of variance to (a) discover "out-rider" values 
in such sets of data, so that they can be charged to "mistakes" 
and deleted so as to prevent an undue influence of these 
values upon the smoothing process, and (b) to formulate 
best-fit smooth functional curves, calculated to minimize 
the sum of squares of the departures of the measured data 
points from the calculated (fitted) curves. 

In the present instance, an additional requirement 
upon the curve-fitting analysis is that it provides the best 
possible cujrve forms for use in extending the domain beyond 
that of the measurements. This is in no way to be construed 
as being equivalent to making measurements in the domain 
extension: it is only an approximation technique that links 
the estimated values in the non-measured region as firmly as 
possible to the measured values. The need for suitable 
laboratory measurements on HCl^g at temperatures below 
273°K is thus emphasized, not met, by the present effort. 

A. Form of Curves 

Particularly in the event that extrapolations are 
required, it is basic that statistical curve fitting be 
guided by physical theory. The form of the dependency of 
vapor pressure on temperature for pure vapors in equilibrium 
with their liquid or solid phases is indicated by the Clausius- 
Clapeyron relation. 
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( 1 ) 


dT R 

V T 

Although L, the latent heat associated with the phase transi- 
tion is quite constant, it can be shown that, if the vapor 
behaves approximately as an ideal gas, L depends upon T 


L = L + £T, 


£ = = c - c 


pv 


( 2 ) 


where c- is the specific heat of the liquid or solid phase 
and Cp^ is the specific heat at constant pressure of the 
vapor. Hence the definite integral of (1) over a temperature 
interval gives an equation of the form 


In p] = A + Bt“^ + C In T 

y 


(3) 


for solutions of constant molality, y. 

The dependency of the equilibrium vapor pressure upon 
concentration is indicated for constant temperature by 
Raoult's Law to have the form 


In p]^ = D + Ey 


(4) 


To provide for departure from ideal conditions, the right 
side of this equation may be extended with terms in powers or 
other functions of y(i.e.. In y, exp y, etc.). In the absence 
of physical reasons for specific choices, the statistical 
procedure of multivariate analysis permits selection of those 
functions that most directly reduce the sum of squares of 
the departures. 
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Obviously the vapor pressures in the present case are 
dependent upon both temperature and molality, hence they may 
most directly be expressed as a surface in (p,T,y) space by 
means of the polynomial 

In p = (A + BT“^ + C In T) (D + + . . . ) (5) 

In the present instance, however, the need to extrapolate 
values to negative Celsius temperatures is basic, and can be 
done most readily on the strength of the Clausius-Clapeyron 
relation (1) and (3). 

B. Procedure 

The water vapor pressure table was considered first 
using equation (3) for the basic form -relating water vapor 
pressure and temperature at each HCl^g concentration listed. 

The least squares curve fitting routine thus provided a set 
of curves: 

which is most readily represented by the array of Aq,Bq, Cq 
values : 


00 

®00 

o 

O 

U 

01 

®01 

^01 


^02 ®02 ^02 


• 

• 

• 

• 

• 

• 

Oi 

®0i 

Si 


for Uq = 0, respectively. 
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L/1 


In the initial step of the analysis, the presence of 
outriders became evident. Three criteria were used to remove 
subjectivity from the decision to classify a value as an 
outrider and delete it from the input data; (a) the departure 
of the measured data point from the smooth p vs T curve 
derived by using all data for each y-value; (b) the cimount 
of reduction of the root mean square error of the fitted 
curve obtained by eliminating the suspect data point; and (c) 
the regularity of the variation of the coefficients A, B and 
C across the array. The following three points were classed 
as outriders and deleted as a result;y = 8.0, T = 273°K; 
y= 11.0, T = 273°K; and y= 15.0, T = 313°K. Unfortunately 
two of these points are at the lower T-limit of the measured 
data, i.e., at the threshold of the region to be extrapolated. 
Results at molalities of 2.0 and 9.0 were relatively poor, but 
no single point qualified as an outrider according to the 
criteria cited above. 

Table 4JT gives the numerical values for the y^ and the 
coefficients, and Figures 41, 42 and 43 show graphically the 
variations of the coefficients as functions of the molality. 

There is no clear physical criterion for the behavior 
of the coefficients A, B and C as functions of y . On the 
other hand, it appears that a non-smooth or discontinuous 
variation of these numbers with molality must be physically 
unrealistic. It is therefore appropriate to smooth A(y), 

B(y), and C(y) using the least squares curve-fitting procedures 
described above. Inasmuch as C (y ) exhibits the least variability 
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Table 4.1: Numerical Values of the Coefficient 

A^, Bq, Cq for Molalities 0 V £ 15.88 


0 

54.413 

-4.9977 

-6789.2 

0.05 

52.489 

-4.7118 

-6701.8 

0.10 

52-729 

-4.7478 

-6712.7 

0.20 

52.262 

-4.6793 

-6690.8 

0.50 

52.228 

-4.6749 

-6691.1 

1.0 

46.635 

-3.8416 

-6445.9 

2.0 

35.498 

-2.1624 

-5989.8 

3.0 

43.796 

-3.4242 

-6338.9 

4.0 

42.654 

-3.2588 

-6301.2 

5.0 

40.248 

-2-9020 

-6215.7 

6,0 

44.018 

-3.4650 

-6411.9 

7.0 

40-240 

-2,9032 

-6272.3 

8.0 

62.165 

-6,1702 

-7293.4 

9.0 

46.358 

-3.8353 

-6582.3 

10-0 

49.957 

-4.3601 

-6804.7 

11.0 

86.551 

-9.8109 

-8496.4 

12.0 

36.227 

-2.3107 

-6275.7 

13.0 

67.505 

-6.9558 

-7751.6 

14.0 

27.670 

-1.0224 

-6998.7 

15.0 

31.200 

-1.5449 

-6207.0 

15.88 

-55.010 

11.424 

-2569.8 
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(Figure 41) from a smooth curve, it was smoothed first leading 
to the revised function C^(y): 

Cj^(y) = -6753.4 + 363. 79y - 66.826y^ + 

' 3 -3 

2.659iy’ + 0.56202 x 10 ^ exp(y). 

New values Aj^(y) and Bj^(y) were then computed using Cj^(y) 
in conjunction with the measured data points, outriders 
deleted, and Bj^(y) was smoothed, giving a revised function 
B2 (y ) : 

B 2 (y) = - 4.8693 + 1.2115y - 0.21394y^ + 

0.86707 X 10“^y^ + 0.19812 x lO"^ exp(y). 

Again new values A 2 (y) were computed using C^(y) and B 2 (y) 
in the regression calculation. Smoothing of these gave the 
revised function A^ (y) : 

A 2 (y) = 53.560 - 8.1475y + 1.4352y^ - 

0.058161y^ - 0.13170 x 10“^ exp(y). 
and the smoothed functions are used in 


In Pjj " A^(y) + B 2 (y)ln T + Cj^(y)T 


-1 


( 6 ) 


to represent the entire water vapor pressure surface in (p,y,T) 
space. Equation (6) predicts the p (y,T) values with 

H2O 

an overall average error of 0.5 percent with respect to the 
tabulated data points. The greatest error of 3.22 percent 
occurs at y= 9.0, T = 273°K, a point at which the outrider 
criteria were nearly, but not quite, fulfilled. 

For ® similar procedure was followed, except 

that 
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MOLALITY 

Fig. 4.1; Coefficienta Aq vs y and the Smoothed 
Curve A^ (y ) . 
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( 7 ) 


In p = A» (y) + B' (y)T~^ + C (y)T 
HCl 

was used Instead of (3) . In turn, the coefficients are related 
somewhat differently to y (Figures 4.4, 4.5, .4.6). 

In this case Aj^(y) was computed first by smoothing the 
data of Figure 44 obtaining 

AJ^(y) = 30.542 - 1.3279y + 0.03242y^ + 

1.6501 In y 

then computing B^<y), C^(y) and smoothing the former to get 

B^(y) = -14694 + 987. 61y - 24.308y^ + 55.001 In y 

and again recomputing to get C^(y) and smoothing to get 

C^(y) = 910.16 X 10^ - 158.^0 x 10^ + 2045. 9y^ 

+ 6331.4 In y + 23439y In y. 

The final form for the equilibrium HCl vapor pressure is 

In = A^(y) + B^(y)T"^ + C^(y)T"^ (8) 

The overall average error given by this formulation as com- 
pared to the Fritz and Fuget (1956) values is about 1 per cent. 
The worst predictions are found at the low molalities and 
reach up to 4.6 percent error. This is not as good as it 
probably should be, and the study is continuing to determine 
whether a better choice of functional parameters can be made. 
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MOLALITY 

6: Coefficients C vs i 

Curve C^ (y ) . 
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